% To calculate casuallity analysis of up- and downgrade std
% transition prob (cohort method) of G1000k annual data
% Reference: 
% 1. de Meulemeester, J.-L., and D. Rochat. 1995. A causality analysis of the link between higher education and economic development. Economics of Education Review 14 (4):351-361.
% 2. http://www.scholarpedia.org/article/Granger_causality
% 3. http://www.spatial-econometrics.com/


clear all;
clc;
load duration1quarters_TP06QRJfinR5G1000kMR7.mat NUMBER_OF_YEAR Quarter_LIST_LENGTH RATE_LIST_LENGTH durationCount
GDP = xlsread('D:\My Documents\CIRANO\Stata\GDP_Canada_quarter.xlsx','c16:c63');

TM = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
xx = zeros (RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        
        x(:, :) = durationCount(i,j,:,:);
        xx = xx + x;  
        
        for d = 1: RATE_LIST_LENGTH-1            
            TM (4*(i-1)+j, d, :) = x(d,:)./ sum(x(d,:),2);            
        end
                
    end
end

TimeLength = size(TM,1);
RateSize = size(TM,2);

% up- and downgrade prob involving a migration of at most two buckets
% downgrade porb: d(k,t) = p(k,k+1,t) + p(k, k+2,t)
% upgrade porb:   u(k,t) = p(k,k-1,t) + p(k, k-2,t)
% except for the extreme categories, where the number of buckets is one)

DownProb = zeros(TimeLength, RateSize);
for yr = 1: TimeLength
    for rate = 1: RateSize
        if rate == RateSize
            DownProb(yr,rate) = TM(yr, rate, rate+1);
        else
            DownProb(yr,rate) = TM(yr, rate, rate+1)+ TM(yr, rate, rate+2);
        end
    end
end

UpProb = zeros(TimeLength, RateSize-1);
for yr = 1: TimeLength
    for rate = 2: RateSize
        if rate == 2
            UpProb(yr,rate) = TM(yr, rate, rate-1);
        else
            UpProb(yr,rate) = TM(yr, rate, rate-1)+ TM(yr, rate, rate-2);
        end
    end
end
 
% Test seasonality of up- and downgrade probability
% regress downProb on seasonaly dummy : DownProb = c+ a1*d1 + a2*d2+
% a3*d3+a4*d4+ miu*t + Y
% Then test a1 = a2 = a3 = a4

d1 = repmat([1 0 0 0]',[12 1]);
d2 = repmat([0 1 0 0]',[12 1]);
d3 = repmat([0 0 1 0]',[12 1]);
d4 = repmat([0 0 0 1]',[12 1]);

t = [1:1:48]';
b = regress(DownProb(:,1),[ones(48,1), d1, d2, d3,d4, t]);


% Test for the presence of unit root (non-stationarity)



% Well-known "augmented" Dickey-Fuller test
for i = 1: RateSize
[h_down(i),~,stat_down(i),cValue_down(i)] = adftest(DownProb(:,i),'model','TS','lags',0);
if i~=1
[h_up(i),~,stat_up(i),cValue_up(i)] = adftest(UpProb(:,i),'model','TS','lags',0); % The first coloumn is 0, should ignore.
end
end

h_GDP = adftest(GDP,'model','TS','lags',1);

% If stationary in level, then apply Granger-causality testing procedure
% The above result shows h=1, reject null of unit root. So the time series
% are stationary level.
[F,c_v] = granger_cause(DownProb(:,7),GDP,0.95,1);
[F,c_v] = granger_cause(GDP, DownProb(:,7),0.95,2);

% If non-stationary in level, test cointegration. In the case of
% conintegration, we test for causality within framework of an error
% correction representation.



